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Abstract 

The methodology used to implement structural 
sensitivity calculations into a major, general- 
purpose finite-element analysis system (SPAR) is 
described. This implementation includes a 
generalized method for specifying element cross- 
sectional dimensions as design variables that can 
be used in analytically calculating derivatives of 
output quantities from static stress, vibration, 
and buckling analyses for both membrane and bend- 
ing elements. Limited sample results for static 
displacements and stresses are presented to indi- 
cate the advantages of analytically calculating 
response derivatives compared to finite difference 
methods. Continuing developments to implement 
these procedures into an enhanced version of SPAR 
are also discussed. 


Introduction 

General-purpose finite-element structural 
analysis programs are widely used in a variety of 
design applications to calculate response 
quantities such as displacements, stresses, 
vibration frequencies, buckling loads, and mode 
shapes to assess the integrity of a proposed 
structure. 1>2,3 Until recently the dominant 
objective of these finite element structural 
analysis programs has been limited to the accurate 
prediction of such structural behavior for a, given 
design: the designer still relies on manual 

evaluation of design modifications based on 
engineering judgment. However, in the structural 
design context, the objective of structural 
analysis should be broadened to include 
calculating estimates of the sensitivity of struc- 
tural response to changes in the design to provide 
a formal approach to guide design modifications. 
This recommendation is based on the fact that 
structural sensitivities can be calculated in a 
straight-forward manner using finite difference 
methods or by using a more accurate and efficient 
approach involving analytical gradients or deriva- 
tives. The established analytical methods 
and user needs indicate the feasibility and 
desirability of including the additional 
capability to calculate analytical structural 
derivatives as a standard option in general- 
purpose, finite-element analysis programs. This 
sensitivity information has several important uses 
such as: (1) a structural designer-computer 

interaction to determine sensitivity of a design 
to changes in structural properties; (2) adjusting 
properties of a finite element model to obtain 
correlation between analysis and test; (3) 
approximating structural analysis using first 
order Taylor series expansion; and (4) calculating 
constraint gradients needed in optimization 
algorithms. Analytical derivatives have been 
incorporated into structural optimization systems 


but such systems often lack generality and/or have 
limited analysis capability. Self-contained 
additions have been made to the computer code of 
the large-scale, general-purpose finite element 
analysis system SPAr3, to include sensitivity 
calculations for static stress and flutter 
analyses as described in Ref. 11. In Ref. 11, a 
linear relationship is specified between design 
variables and the parameters describing the 
structural model. The present paper describes an 
alternate approach for implementing techniques 
which analytically calculate structural response 
derivatives into SPAR. This approach uses 
innovative, specially-constructed sequences of 
input instructions to perform the desired 
calculations with existing analytical procedures 
in SPAR. The unique features of SPAR which 
facilitate this approach are discussed along with 
the generalized method used to allow a nonlinear 
relation between design variables and structural 
model parameters. The equations for the 
structural response derivatives and methods of 
solution are well documented in the literature but 
are summarized herein as a basis for describing 
their implementation into SPAR. ^“9 Limited 
sample calculations are presented for static 
displacements and stresses to indicate the 
advantages of analytically calculating derivatives 
compared to finite difference methods. Also, 
on-going and planned developments for analytical 
derivative calculation are discussed emphasizing 
the anticipated benefits, of implementing these 
procedures in an enhanced version of SPAR. 


SPAR Characteristics 

The capability to calculate analytical deriv- 
atives is implemented primarily by making innova- 
tive use of existing features of SPAR with 
required supplemental computations performed by 
user written subroutines in separate programs. 

The SPAR system is used in a "black-box" mode in 
that required computations are invoked with 
standard SPAR input commands, hence, no internal 
coding changes are needed. 


Modular Organization 

The organization of the SPAR analysis system 
is shown schematically in figure 1. This system 
is composed of a group of individual modules 
called "processors" which are used in a logical 
sequence to perform a desired analysis. Each pro- 
cessor is designed for a limited, yet distinct and 
complete function and is referred to by a unique 
name as shown at the top of figure 1. The func- 
tions of each of the SPAR processors are given in 
table 1. For conventional analysis, the pro- 
cessors TAB through KG read user input, form 
element matrices, and assemble element matrices 
into system matrices which represent the overall 
stiffness and mass of the structure. Equations 
for static stress analysis are solved using 
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processors INV through PSF. The EIG processor is 
used for eigensolutions and AUS and DCU provide 
general input and matrix arithmetic capability. 

A high degree of modularity is provided in 
SPAR since all data communication between 
processors is handled through a data base complex 
which contains the information generated by or 
used by other processors during a computer run. A 
set of data handling utilities transfers data 
between the processors in central memory of the 
computer and the data base complex on auxiliary 
storage. ^2 These utilities may be used in 
auxiliary programs to communicate data to the SPAR 
processors. 


Flexibility of User Input 

User input to the system is shown in the 
upper left portion of figure 1 as executive 
control commands ([XQT followed by the name of the 
processor to be executed) and related processor 
input data located sequentially on the input 
file. A high degree of user flexibility is 
provided since processors can be called for 
execution in any logical sequence. The basis for 
the calculation of structural response derivatives 
described herein is the development of input 
sequence^ to form and solve the structural 
sensitivity equations. 


Coupling SPAR With Auxiliary Proqrams 

The modular organization and flexible user 
input facilitate the use of SPAR analytical 
procedures for "nonstandard" purposes. Such 
usage, as in the implementation of structural 
sensitivity calculations, often requires coupling 
of SPAR to auxiliary programs as illustrated in 
figure 2. Sequences of SPAR input commands for 
invoking processors required to perform desired 
calculations can be prepared manually or generated 
by auxiliary programs. Calculations outside the 
capabilities of the SPAR processors can be 
performed in external programs which use the SPAR 
data handling utilities to communicate with the 
data base complex. Computational sequences 
involving combined and/or repeated executions of 
both SPAR and external programs can be specified 
using the computer operating system control 
language as shown at the top of figure 2. 

Specific programs and procedures used to implement 
the structural response sensitivity calculations 
are presented in subsequent sections of this 
paper. 


Structural Mod i f i cation 


The first step in calculating response 
derivatives for an existing structure is to 
identify what modifications to the finite element 
model are of interest. These modifications are 
defined in terms of design variables such as the 
area or thickness of structural members. The 
first part of this section describes the 
definition of design variables and their 
relationship to the structural model which must be 
established. The second part of this section 
discusses a procedure to produce derivatives of 
the overall stiffness and mass matrices with 
respect to the design variables needed in the 
calculation of structural response derivatives. 


Desigrt Variables Definition 

The user specifies structural modifications 
using design variables, Vi , which in turn must 
be related to parameters which define the 
structural model. Typical parameters which define 
the structural model are section properties or 
mass properties of finite elements in the 
structure. Herein, a structural definition 
parameter, Pj, will be defined to be any 
parameter which has a linear relationship to the 
stiffness and mass matrices of individual finite 
elements in the structural model. Sometimes 
design variables, v, , are defined to be 
identical or in a one-to-one relationship with the 
structural definition parameters, Pj. For 
example, a design variable may be tne area of an 
axial member and the areas of all axial finite 
elements in the structure are defined to be equal 
to that design variable. The relationship between 
Vi and Pj can be more complicated, for 
example, when design variable linking is used or 
bending elements are considered. 

A method for handling a general, nonlinear 
relationship between the structural design 
variables and structural definition parameters 
is presented in Refs. 13 and 14. This generalized 
method of design variable definition requires only 
that derivatives of Pj with respect 
to Vi exist since they are used in the 
calculation of the stiffness and mass matrix 
derivatives. An example of the nonlinearity 
occurs if the moment of inertia per unit width of 
a plate, = t^/12, is a structural parameter 
and if the plate thickness, t, is used as a design 
variable, the derivative 3Pj/8vi is given by 
t 2 / 4 . In another case where the thickness of a 
membrane element, t, is the structural parameter 
and a reciprocal design variable, 1/t is used, the 
derivative, aPj/av^ is given by -t^. Also, 
the derivative, 3Pj/8vi includes design 
variable linking functions when the coefficients 
of such functions are used as design variables. 


Stiffness and Mass Matrix Derivatives 


The derivatives of the stiffness, [K], and 
mass, [M], matrices for the structural model with 
respect to the design variables, v-j , are needed 
in the calculation of structural response 
derivatives. These derivatives are expressed by 
chain differentiation as 


and 
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where NP is the number of structural definition 
parameters. 12 since [K] and [M] are 
linear functions of Pj, the matrix portions of 
the products in equations (1) and (2), 

3[K]/3Pj and 3[M]/3Pj, are constant. 
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These terms can be calculated and saved for 
subsequent use In an Iterative optimization 
procedure. The scalars 3Pj/8yi are not 
necessarily constant and must be calculated for 
particular values of Vi. Formulation by the 
user of equations 1 and 2 for a particular 
application permits consideration of a wide 
variety of structural design problems. 

An illustrative example is a beam whose 
structural definition parameters are the cross- 
sectional area A and moments of inertia 
I2, and do* ^^6 stiffness matrix derivative 
is given by the equation shown at the top of 
figure 3. The structural definition parameters 
are also given for a channel section as nonlinear 
functions of the cross-sectional dimensions. The 
derivatives of the structural parameters at a 
particular value of the design variable are 
obtained by another chain differentiation 
expression as shown for Ij at the bottom of 
figure 3. In general any of the individual 
dimensions may be specified as a design variable. 


Implementation Procedure 

The implementation procedure used to 
calculate stiffness and mass matrix derivatives is 
identical and is shown schematically in figure 4 
for the stiffness matrix. Special sequences of 
SPAR input commands are executed to form the 
desired derivatives as indicated at the bottom of 
the figure. To automate this procedure, auxiliary 
computer programs are used to generate the SPAR 
input sequences. This input can become lengthy 
for large numbers of design variables, v^-, and 
related structural definition parameters, Pj. 

The initial step in this process is to form 
the terms, 3[K]/3Pj and 3[M]/3Pj, as shown 
at the left of figure 4. These terms are formed 
by generating stiffness and mass matrices for unit 
values of the various structural parameters, 

Pj. An example sequence of required SPAR 
input is shown in table 2 for the stiffness 
parameter, Ij, of some beam elements. Repeated 
sets of input, all similar to the example given in 
table 2, are required; each set corresponds to a 
structural parameter for a group of elements that 
undergo the same change in value of structural 
parameter when a design variable is changed. 

In TAB (figure 1) only the unit values of 
structural parameters need be input since tables 
of joint locations, material properties, etc. can 
be used from the data generated during the initial 
definition of the model. The option to specify 
beam properties in terms of stiffness parameters 
is used for setting Ij to unity and all other 
stiffness parameters to zero as shown in table 2. 
The group of elements corresponding to a 
particular structural parameter are input to the 
ELD processor. Processors E and EKS generate the 
unit elemental stiffness matrices, and processors 
K, or M assemble the desired derivative matrices 
3[K]/3Pj and 3[M]/3Pj. These matrices are 
saved in the SPAR data base complex for use in 
subsequent calculations. 

The next step, shown at the right of figure 
4, is to evaluate the scalar terms 3Pj/3Vi 
and to combine them with 3[K]/3Pj and 
3[M]/3Pj as shown in equations 1 and 2 to form 


the stiffness and mass matrix derivatives. The 
terms 3Pj/3vi are dependent on the design 
variable definitions used for a particular 
application. In general, these terms represent 
the evaluation of nonlinear expressions at a 
particular value of the design variable, v^ , and 
would require reevaluation during each iteration 
of a design process. These nonlinear expressions 
are evaluated in user supplied subroutines in an 
auxiliary program and the particular numerical 
values are included in input commands to the AUS 
processor in SPAR where the products and 
summations indicated in equations 1 and 2 are 
performed. An example sequence of SPAR input 
commands to calculate 3[K]/3vi and 3[M]/3v^ 
for a set of beam elements is given in table 3. 
These derivatives are stored in the SPAR data base 
complex under unique names for retrieval and are 
used in subsequent calculations of structural 
response derivatives. 

This procedure for generating the mass and 
stiffness matrix derivatives is satisfactory for 
cases in which many elements are governed by the 
same design variable and a few matrices 
3[K]/3Pj are formed for many elements at a time 
and are multiplied by a single factor 
3Pj/3Vi'. However, this procedure is 
inefficient when, in the most extreme case, 
3[K]/3Pj must be formed for all elements one at 
a time and each subsequently multiplied by a 
different value of 3Pj/3vi. Capabilities of 
a recent, enhanced version of SPAR that offer an 
alternate, improved method for generating the 
stiffness and mass matrix derivatives which 
overcomes this difficulty are discussed in the 
section entitled "Continuing Developments." 


Static Analysis 

Methodology for calculating derivatives of 
displacements and stresses from a static analysis 
are presented in this section. The equations are 
obtained by taking analytical derivatives, with 
respect to specified design variables, of the 
governing structural response equations. 


Derivative Equations 

The matrix equation for static displacements, 
{x}, is 

[K]{x} = {F} (3) 

where [K] is the overall stiffness matrix and 
{F} is a vector of applied loads. Equation 3 is 
differentiated with respect to a design variable 
v^ giving 


CK} 


li x] _ _ i£K]|x} + n 


3Vi 


3Vi 


3Vi 


( 4 ) 


The applied load {F} is usually independent of 
the design variables and 3{F}/3vi is taken to 
be zero herein. The remaining quantity on the 
right hand side of equation 4 is obtained by 
multiplying 3[K]/3Vi from equation 1 by the 
displacement vector {x}. Equation 4 can be 
solved by the same solution algorithm used for 


3 



equation 3, taking advantage of the fact that [K] 
is available in factored form from the solution of 
equation 3. The procedure described above gene- 
rates derivatives for all displacement components 
but requires the solution of equation 4 for each 
design variable. 


A procedure with increased efficiency has 
been implemented for cases where derivatives of a 
selected subset of displacements are 
required. Equation 4 is rewritten as 


ilxL = 


3Vi 


CK]'^{- 





and advantage is taken of the fact that selected 
rows or columns of [K]"l are formed by calcu- 
lating displacements, {q}, due to unit loads, (U}, 
applied at the selected degrees of freedom. A 
significant computational saving results when the 
subset of displacement derivatives is small com- 
pared to the number of design variables, since the 
number of required solutions to the equation 
[K]{q} = {U} is equal to the number of displace- 
ments of interest. The desired derivatives are 
then obtained from 



Element stresses {o} are related to the 
joint displacements by the equation 

{«} = CS]{x} (7) 

and the general expression for stress derivatives 
is 


ii£i 

3V. 
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( 8 ) 


For membrane elements, the matrices of the stress- 
displacement relationships, [S], are independent 
of the design variables and 3[S]/3v-j is zero. 
Therefore, stress derivatives can be calculated by 
simply repeating the stress calculation procedure, 
equation 7, with the displacement derivatives 
3{x}/3v-j replacing the displacements {x}. 

For bending elements, the stress-displacement 
relationships are dependent on the element cross 
section geometry and this dependence must be 
included in the stress derivative calculations. 

As an example, the stresses for a beam are given 
by 


F, M, 

°p = 7T + T7^p2 - 1^ ^pl 

where p refers to a specified point of the beam 
cross section. F3 is the longitudinal force 
acting at the centroid of the cross sectional 
area. A; and M3, M2, 13, and I2 are 
applied moments and area moments of inertia about 


the principal axes of the beam cross section. The 
distance of the specified point from the principal 
axes are given by yp2 and yp3. 

Differentiating equation 9 with respect to a 


design 

variable, v-j , 

gives 
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The terms in the square brackets, [ ], can be 
calculated using geometric properties of the beam 
cross section and are multiplied by F3, M3, 
and M2 from a static analysis of the original 
structure. Evaluation of the force and moment 
derivatives 3F3/3vi, 3M3/3vi, and 
3M2/3v-j is performed by repeating the 
procedure for calculating element forces and 
moments with the displacement derivatives 
3 {x}/ 3Vi' replacing the displacements {x| and 
these derivatives are subsequently multiplied by 
the geometric terms as shown in equation 10. 


Implementation Procedure 

A static stress analysis of the original 
structural model is performed as the first step in 
calculating displacement and stress gradients. 

The SPAR processors TAB through PSF are used for 

these calculations. The resulting displacement 

vectors for all applied load cases are multiplied, 

using processor AUS, by all previously generated 

stiffness matrix derivatives to form the set of 

pseudo-load vectors shown on the right side of 

equation 4. These pseudo-load vectors are given 

names corresponding to a set of applied forces 

with unique load set numbers. Two options exist 

for solving for the displacement gradients. In 

option I, for use when the number of desired 

displacement gradient components is larger than 

the number of design variables, the SSOL processor 

is used to solve equation 4 directly. In option 

II, for use when the number of design variables is 

larger than the number of desired displacement 

gradient components, the SSOL processor is used to 

calculate displacements (q} due to unit loads at 

the selected degrees of freedom. The displacement 

derivatives shown in equation 6 are then formed in • 

AUS by multiplying {q}'’’ by the pseudo-load 

vectors. 

Required stress derivative information for 
membrane elements is generated using the 6SF 
processor with the displacement gradients 
specified as input in place of the actual 
displacements as indicated in equation 8. For 
beam elements, an auxiliary program which uses 
SPAR data directly as illustrated at the right of 
figure 2 is used to evaluate the stress 
derivatives. This program retrieves previously 
generated forces and moments and their derivatives 
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from the SPAR data base complex and combines them 
with values produced by subroutines which evaluate 
the geometric terms given in equation 10. The 
SPAR data handling utilities are used in the 
auxiliary program to store the resulting beam 
stress derivatives in the SPAR data base complex. 


Vibration and Buckling Analysis 

Methodology for calculating derivatives of 
eigenvalues and eigenvectors is presented in this 
section. The equations for vibration analysis are 
presented and modifications required for buckling 
analysis are indicated. 


Derivative Equations 

The matrix equation for vibration analysis is 
[K-u2m]{4i} = 0 (11) 


The mode and frequency derivatives are both 
obtained by the solution of equation 16. This 
method requires the formation of the matrix on the 
left side of equation 16 which has an additional 
row and column and hence a different topology than 
the original global stiffness and mass matrices. 
Since the SPAR system handles the global matrices 
in a special sparse format, the addition of a row 
and column would require the judicious use of 
fictitious elements with properties selected to 
■give the desired terms in the equations. 

An approximate method for calculating the 
eigenvector derivatives, described in Ref. 4, is 
available in SPAR as part of the SM processor 
which is used to modify the properties of 
structural models to correlate analytical and test 
results as discussed in Ref. 16. In this method, 
the eigenvector derivatives are approximated as a 
linear combination of a subset of eigenvectors of 
the original structure. The derivative of the jth 
eigenvector is expressed as 


where and {<(i} are sets of eigenvalues and 
eigenvectors corresponding to the natural 
frequencies and mode shapes of the structural 
model. The normalization of the modes with 
respect to the mass matrix is given by 

= [I] (12) 

Differentiating equation 11 with respect to the 
design variable v-j gives 


[K.„2M]iilL = lit 

•'3V^ 3V. 
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and the corresponding differentiation of equation 
12 gives 


(14) 


■*av. 


If only derivatives of the frequency are 
required, equation 13 is premultiplied by 
and when it is recognized that 
equals 0.0 and {<ti}*3(o2/3vi [M]{(|i equals 
3u)2/3vi the equation becomes 


%■ - 




(15) 


Several methods exist to solve for the 
derivatives of the natural vibration modes 
3()i/3Vi. ^>7,8 jhe method presented in Ref. 8 

combines equations 13 and 14 in the form 
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(17) 


The contribution of the kth eigenvector is 
determined by substituting the expression for 
3{i(>j}/3vi into equation 13 and premultiplying 
by UkT to give 


(18) 


and by a similar substitution into equation 14 to 
give 

aijk = ^ for k=j (19) 


The disadvantage of this method is the uncertainty 
of how many eigenvectors are required to give an 
adequate representation of the eigenvector 
derivatives. 


The method for calculating eigenvector 
derivatives that is described in Ref. 7 has been 
implemented using sequences of SPAR input commands 
for performing the required calculations. In this 
method, the right side of equation 13 is evaluated 
and treated as a pseudo-load vector. A direct 
solution of equation 13 is not possible 
since [K-uj2m] is a singular matrix. As 
discussed in Ref. 7, if the value of one component 
of 3{ii)}/3v-j is fixed, the remaining components 
can then be calculated yielding a solution {P} to 
equation 13. Since the eigenvector {(t>} is the 
homogeneous solution of equation 13, the following 
expression is also a solution 

{P} (20) 
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The value of the multiplier, C, for which equation 
20 also satisfies equation 14 is given by 

C = (21) 

The equations for buckling response derivatives 
are similar, with the geometric stiffness matrix 
replacing the mass matrix. 


Implementation Procedure 

The procedure for calculating vibration 
response derivatives, described in Ref. 7, was 
implemented by manually preparing the required 
SPAR input commands to demonstrate the method. 
Computer programs for automatic generation of this 
input similar to those for static analysis 
derivatives have not been developed. Use is made 
of the static analysis procedures, discussed in 
previous sections, to generate the stiffness and 
mass matrices, [K] and [M], and their 
derivatives, 3[K]/3vi and 3[M]/3Vi. The EIG 
processor in SPAR is used to calculate 
eigenvalues, and eigenvectors, {(ti}, for the 
original structure. These terms are combined, 
using the AUS processor, to form the frequency 
derivatives shown in equation 15, The right side 
of equation 13 and the matrix [K-w^m] in 
sparse format is also formed using the AUS 
processor. A solution, {P}, to equation 13 is 
calculated using the INV and SSOL processors with 
one component of {P} set to zero using a 
fictitious boundary constraint. The component 
corresponding to the largest component 
of (iij} is selected following the 
recommendation of Ref. 7 based on numerical 
accuracy considerations. Finally, the required 
vibration mode derivatives are calculated using 
the AUS processor to perform the matrix algebra 
for evaluating equations 21 and 20, respectively. 
The entire procedure must be repeated for each 
eigenvector derivative. 

, In the calculation of buckling response 
derivatives, similar steps are performed with the 
KG processor used to calculate the differential 
stiffness matrix [Kg]. The stress derivatives, 
obtained from a previously described procedure are 
input to the KG processor to form the 
sTKgl/svi matrix that is required in the 
buckling derivative equations. 


Sample Calculations 

The analytical gradient capabilities that are 
implemented in SPAR were used to calculate 
derivatives for comparison with similar results 
from finite difference methods. The use of 
analytical gradients provides improved 
computational efficiency and eliminates the 
uncertainty of numerical accuracy associated with 
selection of a finite difference perturbation to 
give desired convergence. 


Model Description 

The finite element structural model of a 
stiffened cylinder containing a rectangular 


cutout, as shown in figure 5, is used to 
illustrate the calculation of static displacement 
and stress derivatives. ] Results are calculated 
for three different levels of model refinement to 
assess the effect of number of degrees-of-freedom 
on the required computational time. The 337 
degree-of-freedom cylinder model, referred to as 
Model 1, is stiffened by five rings equally spaced 
along its length and sixteen stringers equally 
spaced around its ci rcumference as illustrated in 
the figure. The rings and stringers are modeled 
by beam and rod elements respectively. 

Rectangular panels between rings and stringers are 
modeled by membrane plate elements. The cutout is 
located between the second and fourth rings and 
encompasses a 90° segment of the cylinder's 
circumference. The translational degrees of 
freedom are constrained at all points on one end 
of the cylinder. The level of modeling is refined 
for Models 2 and 3 by dividing the skin panels 
with additional rings and stringers to give the 
total number shown on figure 5. Two sets of 
concentrated forces are applied at the unsupported 
end of the cylinder as shown in the figure. 

Three design variables are considered: (1) 

an area which governs all stringer rod elements; 
(2) a thickness for all membrane plates; (3) a 
scale factor on the dimensions of a specified 
chahnel cross-section used for all beams in the 
rings. Hence, each design variable affects many 
elements in the same manner and results in a 
simplified calculation of the derivatives, 
3[K]/3Pj, 


N umerical Results 

The computation times required for each of 
the major steps in calculating stress and 
displacement derivatives for the two load cases 
are shown in figure 6 for the three levels of 
model refinement. Performing a static analysis 
with processors INV and SSOL is the step that 
requires the largest amount of time. This time is 
used primarily by INV to factor the stiffness 
matrix. However, only a single static analysis is 
needed rather than multiple analyses for 
perturbations of all the design variables one at a 
time as required by the finite difference method. 
The times shown in figure 6 are CPU (central 
processing unit) seconds on a CDC Cyber 173 
computer when three design variables and two 
loading cases are considered. 

The relative computer time (finite difference 
method divided by analytical gradient method) is 
shown in figure 7. The results for the case of 
three design variables are extrapolated to show 
the effect of a larger number of design 
variables. This extrapolation is based on 
appropriate scaling of the times shown in figure 6 
to reflect the number of times each step in the 
process must be repeated for additional design 
variables. The relative times and hence benefits 
increase with the number of degrees of freedom. 
Benefits of the analytical gradient method also 
increase with an increase in the number of design 
variables, especially for the most refined model. 
These results were obtained using the procedure 
which calculates displacement derivatives for all 
degrees of freedom and stress derivatives for all 
elements. 
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Use In Optimization Procedures 

The methods used to form and solve the 
analytical derivative equations are organized for 
efficient use in structural optimization 
systems. The major steps in an optimization 
procedure which includes displacement and stress 
constraints are shown schematically in figure 8. 
All steps required to compute response 
sensitivities or gradients are shown in the left 
and center of the figure. The optimization loop 
to update the design variables is indicated at the 
right of the figure. The steps or operations 
shown in the solid boxes are performed by 
executing SPAR with previously described input 
sequences. Auxiliary programs are used to perform 
the operations in the dashed boxes including 
calculation of 3Pj73vi using problem 
dependent subroutines and corresponding 
subroutines needed to calculate the beam stress 
derivatives. The formation of the linear 
derivatives 3[K]/3Pj is not included in the 
optimization loop since it is performed as a 
preprocessor operation and need not be repeated. 


Continuing Developments 

A recently enhanced version of SPAR, called 
EAL (Engineering Analysis Language)^^ permits 
improvements in the implementation of the 
structural response sensitivity calculations. The 
existing procedures for gradient calculation which 
have been discussed herein are being implemented 
in EAL. Current efforts and plans for this work 
are discussed in this section. 

The enhancements in EAL which make the 
improvements in gradient calculation capability 
possible are shown in table 4. Facilities are 
provided for addition of new processors to allow 
all computational procedures which were previously 
handled in auxiliary programs to be integrated 
into the EAL system with all data communication 
through a common data base complex. The addition 
of new processors in the EAL system is similar to 
the approach taken in Ref. 11. 

Sequences of input instructions can be stored 
as data sets in the EAL data base complex and 
retrieved by name for execution. This feature 
provides a more unified implementation of the 
gradient capability than defining the input 
sequences for various steps in the gradient 
calculations as external files to be handled by 
the computer operating system. Also, looping and 
branching statements available in EAL permit 
repetitive execution of a sequence of such steps 
for calculation of vibration mode derivatives for 
each natural frequency or for steps in an 
iterative optimization procedure. 

A processor called LSK is available in EAL to 
assemble a partial system stiffness matrix which 
contains contributions from elemental stiffness 
matrices for a user-selected set of elements. Use 
of this processor simplifies the calculation of 
the stiffness matrix derivatives 3[K]/av-j. The 
terms 3[K]/aPj are formed for all elements at 
once by generating individual stiffness matrices 
with unit values of the structural parameters. 

The LSK processor can be used subsequently to 
scale and assemble these matrices to produce the 
stiffness matrix derivatives a[K]/3v-j as shown 


in equation 1. A table is input to LSK to specify 
which subset of element matrices are to be 
included during an execution and another table is 
input to specify a multiplying factor to be 
applied to the stiffness matrix of each element as 
it is being assembled. The first table can be 
used to simplify the general specification of 
design variable linking with the values in the 
second table corresponding to the aPj/avi 
terms in equation 1. Another approach which has 
proved to be effective is to use the LSK processor 
to retrieve element stiffness matrices needed to 
calculate the stiffness matrix derivatives using 
finite difference methods. The finite difference 
calculation of 3[K]/avi allows any structural 
model definition input parameter to be used as a 
design variable. 

Implementation of procedures using EAL to 
produce gradients of structural response with 
respect to changes in joint locations of the model 
has been investigated. Finite difference methods 
are used to calculate the stiffness and mass 
matrix derivatives, a[K]/av-( and a[M]/avi, 
and the remainder of the steps are the same as 
those where member sizes are designated as design 
variables. The finite difference perturbation of 
joint locations is expedited by use of parametric 
definition of model geometry which is available in 
EAL. Variables, whose values are set at execution 
time, are used in the input instructions along 
with the joint mesh generation capabilities of EAL 
to define the parametric models. For example, all 
joint locations in a wing structure could be 
changed by changing the variable defining sweep 
angle. 

Continued utilization of enhancements 
available in EAL is leading to an integrated, 
large-scale capability to calculate structural 
response sensitivities for a variety of design 
variables. 


Concluding Remarks 

The methodology used to implement structural 
sensitivity calculations into a major, general- 
purpose finite element analysis system (SPAR) was 
described. This implementation includes a 
generalized method for specifying element cross- 
sectional dimensions as design variables that can 
be used in analytically calculating derivatives of 
output quantities from static stress, vibration, 
and buckling analyses for both membrane and 
bending elements. Limited sample results for 
static displacements and stresses were presented 
to indicate the advantages of analytically 
calculating response derivatives compared to 
finite difference methods. Continuing 
developments to implement these procedures into an 
enhanced version of SPAR were also discussed. 
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I mprovements 

Integrated system with all data coimuni- 
cated internally through data base 

Replace use of computer operating system 
control language 


Simplify existing method of calculating 
3[K]/8Vi 

Finite-difference calculation of 
3[K)/9Vi 


Parametric definition of model geometry 


Enhancements 

Facilities for adding new processors 

Input sequences stored in and called 
from data base 

Looping and branching commands in 
input sequences 

LSK processor to assemble stiffness 
matrices of selected elements 


Input variables whose values are 
set at execution time 
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Fig. 4 Implementation procedure used to calculate 
stiffness matrix derivatives. 


Fig. 1 Organization of SPAR analysis system. 



Fig. 2 Coupling of SPAR and auxiliary programs 
for nonstandard calculations. 


LOAD 
CASE 1 

20 000 MEMBER SIZE PROPERTIES 



Fig. 5 Example structural model of stiffened 
cylinder with cutout. 
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Fig. 3 Stiffness matrix derivative expressions 
for channel beam cross-sectional changes. 
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Fig. 6 Computation time for steps required to 
calculate static displacement and stress 
derivatives. 
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Fig. 7 Relative computer time for calculation of 
static derivatives using finite difference 
and analytical gradient methods. 



Fig. 8 Implementation of analytical structural 
response derivatives in a structural 
optimization procedure. 
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